#ifndef __CRAZYFISH_MATRIX__//避免重复调用本头文件，第二次以后就跳到donothing
#define __CRAZYFISH_MATRIX__

#include <iostream>
#include <new>
#include <vector>
#include <exception>
#include <cmath>
#include <limits>
using namespace std;

/// Set the machine epsilon.
const double eps = numeric_limits<double>::epsilon();

class SquareMatrix
{
private:
    double* matrix_elements;
    int dimension;
    void allocMem(int _d);
    void releaseMem();
    void exitAbnormal(int _err_code);
    void copy_data(const double* _f, double* _t, int _d);
    
public:
    SquareMatrix()//缺省构造函数
	{
	    matrix_elements = NULL;
	    dimension = 0;
	}
    SquareMatrix(int _n);
    SquareMatrix(const SquareMatrix& rhs);
    ~SquareMatrix();
    const SquareMatrix& operator=(const SquareMatrix& rhs);
    double& operator()(int _i, int _j);
    vector<double> invMul(const std::vector<double>& rhs);
    
    friend vector<double> operator/(const vector<double>& lhs, const SquareMatrix& rhs);
};

vector<double> operator/(const vector<double>& lhs, const SquareMatrix& rhs);//二元运算符因对称，常常定义在类外，需要friend声明以调用类内private

void exitAbnormal(int error_code);

void exitAbnormal(int error_code)
{
    switch (error_code)
    {
    case 1:
	cerr << "The system is linear dependence." << endl;
	break;
    case 2:
	cerr << "The dimensions of the vector and the matrix are not match." << endl;
	break;
    default:
	break;
    }
    terminate();
};

void SquareMatrix::copy_data(const double* _f, double* _t, int _d)
{
    for (int i = 0; i < _d; i++)
	_t[i] = _f[i];
};

void SquareMatrix::exitAbnormal(int _err_code)
{
    switch (_err_code)
    {
    case 1:
	cerr << "index out of range." << endl;
	releaseMem();
	break;
    case 9:
	cerr << "no enough memory." << endl;
	break;
    default:
	break;
    }
    terminate();
};

void SquareMatrix::releaseMem()
{
    if (matrix_elements != NULL)
	delete [] matrix_elements;
    matrix_elements = NULL;
};

void SquareMatrix::allocMem(int _d)
{
    try
    {
	matrix_elements = new double[_d];
    }
    catch (bad_alloc& no_memory)
    {
	cerr << no_memory.what() << endl;
	exitAbnormal(9);
    }
};

double& SquareMatrix::operator()(int _i, int _j)
{
    if (_i < 0 || _i > dimension - 1 || _j < 0 || _j > dimension - 1)
	exitAbnormal(1);
    return *(matrix_elements + _i * dimension + _j);
};

const SquareMatrix& SquareMatrix::operator=(const SquareMatrix& rhs)
{
    int n = rhs.dimension;
    if (matrix_elements != NULL)
	releaseMem();
    allocMem(n * n);
    copy_data(rhs.matrix_elements, matrix_elements, n * n);
    dimension = n;
    return *this;     
};

SquareMatrix::SquareMatrix(const SquareMatrix& rhs)
{
    int n = rhs.dimension;
    allocMem(n * n);
    copy_data(rhs.matrix_elements, matrix_elements, n * n);
    dimension = n;
};

SquareMatrix::SquareMatrix(int _n)
{
    allocMem(_n * _n);
    dimension = _n;
};

SquareMatrix::~SquareMatrix()
{
    releaseMem();
};

vector<double> operator/(const vector<double>& lhs, const SquareMatrix& rhs)
{
    int n = rhs.dimension;
    vector<double> x(n);
    SquareMatrix A(rhs);
    vector<double> b(lhs);
    if (n != lhs.size())
	exitAbnormal(2);
    /// A gaussian elimination procedure without pivoting.
    for (int k = 0; k < n - 1; k++)
    {
	/// Check the linear dependency of the system.
	if (fabs(A(k, k)) < eps)
	    exitAbnormal(1);
	for (int i = k + 1; i < n; i++)
	{
	    for (int j = k + 1; j < n; j++)
		A(i, j) -= A(k, j) * A(i, k) / A(k, k);
	    b[i] -= b[k] * A(i, k) / A(k, k);
	}
    }
    if (fabs(A(n - 1, n - 1)) < eps)
	exitAbnormal(1);

    /// Backward substitution.
    for (int i = n - 1; i >= 0; i--)
    {
	double sum = 0.0;
	for (int j = i + 1; j < n; j++)
	    sum += A(i, j) * x[j];
	x[i] = (b[i] - sum) / A(i, i);
    }
    return x;
};

vector<double> SquareMatrix::invMul(const vector<double>& rhs)
{
    int n = dimension;
    vector<double> x(n);
    SquareMatrix A(*this);
    vector<double> b(rhs);
    if (n != rhs.size())
	exitAbnormal(2);
    /// A gaussian elimination procedure without pivoting.
    for (int k = 0; k < n - 1; k++)
    {
	/// Check the linear dependency of the system.
	if (std::fabs(A(k, k)) < eps)
	    exitAbnormal(1);
	for (int i = k + 1; i < n; i++)
	{
	    for (int j = k + 1; j < n; j++)
		A(i, j) -= A(k, j) * A(i, k) / A(k, k);
	    b[i] -= b[k] * A(i, k) / A(k, k);
	}
    }
    if (std::fabs(A(n - 1, n - 1)) < eps)
	exitAbnormal(1);

    /// Backward substitution.
    for (int i = n - 1; i >= 0; i--)
    {
	double sum = 0.0;
	for (int j = i + 1; j < n; j++)
	    sum += A(i, j) * x[j];
	x[i] = (b[i] - sum) / A(i, i);
    }
    return x;
};

#else
/// DO NOTHING.
#endif
